Code
library(tidyverse)
library(here)
library(broom)
# Spatial data packages
library(sf)
library(tmap)<p>Raw denim you probably haven't heard of them jean shorts Austin. Nesciunt tofu stumptown aliqua, retro synth master cleanse. Mustache cliche tempor, williamsburg carles vegan helvetica. Reprehenderit butcher retro keffiyeh dreamcatcher synth. Cosby sweater eu banh mi, qui irure terry richardson ex squid. Aliquip placeat salvia cillum iphone. Seitan aliquip quis cardigan american apparel, butcher voluptate nisi qui.</p>
<p>Food truck fixie locavore, accusamus mcsweeney's marfa nulla single-origin coffee squid. Exercitation +1 labore velit, blog sartorial PBR leggings next level wes anderson artisan four loko farm-to-table craft beer twee. Qui photo booth letterpress, commodo enim craft beer mlkshk aliquip jean shorts ullamco ad vinyl cillum PBR. Homo nostrud organic, assumenda labore aesthetic magna delectus mollit.</p>
<p>Etsy mixtape wayfarers, ethical wes anderson tofu before they sold out mcsweeney's organic lomo retro fanny pack lo-fi farm-to-table readymade. Messenger bag gentrify pitchfork tattooed craft beer, iphone skateboard locavore carles etsy salvia banksy hoodie helvetica. DIY synth PBR banksy irony. Leggings gentrify squid 8-bit cred pitchfork.</p>
<p>Trust fund seitan letterpress, keytar raw denim keffiyeh etsy art party before they sold out master cleanse gluten-free squid scenester freegan cosby sweater. Fanny pack portland seitan DIY, art party locavore wolf cliche high life echo park Austin. Cred vinyl keffiyeh DIY salvia PBR, banh mi before they sold out farm-to-table VHS viral locavore cosby sweater.</p>
library(tidyverse)
library(here)
library(broom)
# Spatial data packages
library(sf)
library(tmap)#read in the CA Counties shapefile
ca_counties_raw_sf <- read_sf(here("a2_task1_lecuona_sophia", "data", "ca_counties_copy", "CA_Counties_TIGER2016.shp"))
# clean names and go from m^2 to km^2 for land area
ca_counties_sf <- ca_counties_raw_sf %>%
janitor::clean_names() %>%
mutate(land_km2 = aland / 1e6) %>%
select(county = name, land_km2)
# ca_counties_sf %>% st_crs()
# ca_counties_sf %>% terra::crs()
# ^ done to see it is currently WGS 84 / Pseudo-Mercator#read in the Oil Spill Incident Tracking shapefile
oil_raw_sf <- read_sf(here("a2_task1_lecuona_sophia", "data", "ca_oil", "Oil_Spill_Incident_Tracking_[ds394].geojson"))
# oil_sf <- oil_raw_sf %>%
# select(-LATITUDE, -LONGITUDE)
write_sf(oil_raw_sf, here("a2_task1_lecuona_sophia", "data", "ca_oil", "oil_tracking.gpkg"))
oil_sf <- read_sf(here("a2_task1_lecuona_sophia", "data", "ca_oil", "oil_tracking.gpkg")) %>%
janitor::clean_names() %>%
select(-latitude, -longitude)
# Check the CRS:
# oil_sf %>% st_crs()
# oil_sf%>% terra::crs()
# ^ currently WGS 84 #exploratory plot
plot(oil_sf %>% select(objectid))# transform the projection to match:
oil_projected_sf <- st_transform(oil_sf, st_crs(ca_counties_sf))
# oil_projected_sf %>% st_crs()
# check it! Yay it worked. Both are WGS 84 / Pseudo-Mercatorggplot() +
geom_sf(data = ca_counties_sf) +
geom_sf(data = oil_projected_sf, size = 1, color = "red")# spatial join!
ca_and_oil_sf <- ca_counties_sf %>%
st_join(oil_projected_sf)
# by counts of incidents
ca_oil_counts_sf <- ca_and_oil_sf %>%
group_by(county) %>%
summarize(n_records = sum(!is.na(objectid)))
# # change 'NA' row to be more intuitive:
# ca_oil_counts_sf[58, 1] <- "Pacific Ocean"
head(ca_oil_counts_sf)Simple feature collection with 6 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -13668380 ymin: 4502606 xmax: -13307390 ymax: 4888059
Projected CRS: WGS 84 / Pseudo-Mercator
# A tibble: 6 × 3
county n_records geometry
<chr> <int> <MULTIPOLYGON [m]>
1 Alameda 188 (((-13612247 4538150, -13612347 4538291, -13612423 453844…
2 Alpine 1 (((-13366504 4678946, -13366493 4678954, -13366468 467897…
3 Amador 8 (((-13472698 4647652, -13472698 4647677, -13472698 464775…
4 Butte 11 (((-13565005 4798394, -13564992 4798529, -13565000 479854…
5 Calaveras 7 (((-13428575 4627725, -13428535 4627889, -13428535 462809…
6 Colusa 4 (((-13589905 4781178, -13589881 4781178, -13589804 478117…
# Set the viewing mode to "interactive":
tmap_mode(mode = "view")
tm_shape(ca_counties_sf) +
tm_fill("land_km2", palette = c("#B0E0E6", "#DAA520","#A0522D")) +
tm_shape(oil_projected_sf) +
tm_dots(col= "darkblue") +
tm_layout(title = "Oil Spill Incident Tracking",
title.size = 1.5,
frame = FALSE)